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STATIONARY DISCRETE SHOCK PROFILES FOR SCALAR 
CONSERVATION LAWS WITH A DISCONTINUOUS GALERKIN 

METHOD 


FLORENT RENAC* 

Abstract. We present an analysis of stationary discrete shock profiles for a discontinuous 
Galerkin method approximating scalar nonlinear hyperbolic conservation laws with a convex flux. 
Using the Godunov method for the numerical flux, we characterize the steady state solutions for 
arbitrary approximation orders and show that they are oscillatory only in one mesh cell and are 
parametrized by the shock strength and its relative position in the cell. In the particular case of the 
inviscid Burgers equation, we derive analytical solutions of the numerical scheme and predict their 
oscillations up to fourth-order of accuracy. Moreover, a linear stability analysis shows that these 
profiles may become unstable at points where the Godunov flux is not differentiable. Theoretical 
and numerical investigations show that these results can be extended to other numerical fluxes. In 
particular, shock profiles are found to vanish exponentially fast from the shock position for some class 
of monotone numerical fluxes and the oscillatory and unstable characters of their solutions present 
strong similarities with that of the Godunov method. 
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1. Introduction. Discontinuous Galerkin (DG) methods are high-order finite 
element discretizations which were introduced in the early 1970s for the numerical 
simulation of the first-order hyperbolic neutron transport equation [iniisn]. In recent 
years, these methods have become very popular for the solution of nonlinear convec¬ 
tion dominated flow problems [aiHiiiT]. These methods allow high-order of accuracy 
and locality, which make them well suited to parallel computing, /ip-refinement, hp- 
multigrid, unstructured meshes, application of boundary conditions, etc. 

However, the DG method suffers from spurious oscillations in the vicinity of 
discontinuities that develop in solutions of hyperbolic systems of conservation laws. 
These oscillations are due to the Gibbs phenomenon [9] and may cause the solution to 
become locally nonphysical leading to robustness issues for the computation. Quadra¬ 
ture rules are usually used to compute integrals in the discretization of the equations. 
Properties of the DG method therefore depend on the local structure of the numerical 
solution at faces and within elements of the mesh. The control of such oscillations 
at a reasonable cost while keeping accuracy, robustness and stability is essential for 
the efficiency of the DG method and remains a challenge. Strategies have been pro¬ 
posed such as limiters 0 0, non-oscillatory reconstructions 0 129) . /ip-adaptation 

shock capturing techniques [iziiin], etc. The latter methods aim at adding arti¬ 
ficial viscosity to spread the structure of the discontinuity so that it can be resolved 
at the discrete level. For a DG method with polynomials of degree p > 0 and cells of 
size h, the resolution scale is h/p thus meaning that the method should resolve the 
discontinuity inside the element m- The behavior of the numerical solution inside 
the discretization elements is therefore important for understanding the convergence 
of the DG computations. 

In this work, we focus on the DG method for scalar nonlinear conservation laws 
which present stationary discontinuous solutions. More precisely, we are interested 
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in the behavior of discrete profiles near shocks. Jennings [13] studied the approx¬ 
imation of scalar equations by monotone conservative finite difference schemes and 
proved the existence and stability of traveling discrete shocks. Discrete profiles for 
scalar equations were also studied in miniiin. The analysis was then extended to 
systems in |24l . Bultelle et al. |4] analyzed the linear stability of steady shock 

prohles obtained for general systems of conservation laws discretized with the Go¬ 
dunov method and constructed unstable profiles in the case of the Euler equations for 
gas dynamics. Recently, Lerat [IH] found exact discrete shock solutions for residual- 
based compact schemes m up to seventh-order of accuracy discretizing the inviscid 
Burgers equation. Solutions were derived explicitly and were parametrized by the 
relative position of the shock in the discretization cell. Such analysis may help to 
tune parameters of the numerical method. In the context of DG methods, Cockburn 
and Guzman |0] considered a formally second order approximation of a scalar linear 
hyperbolic equation with discontinuous initial condition. They gave estimates of the 
size of extent of oscillations upstream and downstream of the discontinuity based on 
suitable weights introduced in m- Recently, their work was extended to arbitrary 
approximation order in space and a third-order Runge-Kutta method on non-uniform 
meshes [M] [3S| . The case of ^-singularities in initial condition and source term of a 
scalar hyperbolic equation was then investigated in |33j where superconvergence in 
negative-order norms outside the pollution region was proved and sharp estimates 
over the whole domain were given. 

The objective of this work is the theoretical analysis of stationary discrete shock 
solutions for the DG method discretizing a scalar conservation law with a general con¬ 
vex flux. We will mainly consider the Godunov method to evaluate the numerical flux 
for which explicit solutions can be derived, but we will also focus on other numerical 
fluxes widely used in the context of DG methods. To this end, we first give general 
results about the structure of the stationary discrete profiles. These results are local 
in the sense that they consider profiles that are small perturbations to the exact solu¬ 
tion. In the case of the Godunov flux, we establish the mean value of the solution in 
the cell containing the exact shock position and the complete solution in other cells 
for an arbitrary approximation order. The analysis also predicts exponential decay 
of oscillations of the discrete profile on both sides of the shock for a certain class of 
monotone numerical fluxes. Then, the linear stability of these profiles is investigated 
and eigenvalues are characterized for the Godunov flux. As a result, the shock prohles 
may become unstable at points where the numerical hux is not differentiable. These 
points contain the situation of a shock at interface which is counter-intuitive since 
the exact solution is a piecewise constant function over all cells and is included in the 
function space of the DG method. 

Gonsidering the inviscid Burgers equation we derive exact discrete shock pro¬ 
hles for the DG scheme up to fourth-order of accuracy in the spirit of the work of 
Lerat [18]. The results allow to evaluate quantitatively the structure of the solution 
within elements and predict situations where the numerical solution violates locally 
an entropy condition at the cell level. The theoretical stability analysis predicts the 
occurrence of unstable prohles when the exact shock position is close enough to an 
interface of the mesh. Numerical experiments will suggest that many features of the 
DG method obtained with the Godunov hux still hold for other numerical huxes. In 
particular, for a given space discretization, the oscillating and unstable characters of 
the solution depend on the strength and exact position of the shock. As this latter fea¬ 
ture is generally unknown, it strongly complicates the analysis for enhancing stability 
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and robustness of the DG method. These results may support approaches based on a 
posteriori limitation techniques such as the MOOD method [ills]. Artificial viscosity 
represent another attractive approach providing that the amount of viscosity adapts 
itself to the regularity of the solution. Theoretical results in this work will suggest the 
addition of artificial viscosity to the highest modes in the DG space spanned by hier¬ 
archical basis functions. As an illustration, we apply the spectral vanishing viscosity 
method jH] with a selective filter in the Legendre basis as proposed in p5] . 

The paper is organized as follows. Section [^ presents the model problem and 
the numerical approach for the discretization. In section we consider a general 
convex flux and analyze discrete shock profiles for the DG method. These profiles are 
explicitly constructed in the case of the inviscid Burgers equation in section|4 A linear 
stability analysis is performed in the vicinity of these solutions in section ^ These 
results are assessed by several numerical experiments in section ^ and a first attempt 
for stabilizing the numerical scheme is proposed in section [673 Finally, concluding 
remarks about this work are given in section 

2. Model problem and discretization. 

2.1. Nonlinear scalar equation. The discussion in this paper focuses on the 
discretization of scalar nonlinear hyperbolic equations in one space dimension with a 
DG method. Let D = K be the space domain and consider the following problem 


dtu + d^f{u) 
u{x, 0) 


0 , 

Uo{x), 


in D X (0, oo), 
in D. 


(2.1a) 

(2.1b) 


The physical flux / in C^(yia) is assumed to be coercive and strictly convex over 
the set of admissible states Da C M: lima_>±oo f{u) = -|-c» and f"{u) > 0 for all u in 
Da. We are particularly interested in steady shock solutions to (13. Such solutions 
consist in stationary discontinuities between two states 


u{x) := lim u{x,t) = 

t—^OO 


ul if X < Xc, 

Ur if X > Xc, 


( 2 . 2 ) 


where Xc 
relation 


denotes the shock position and the states satisfy the Rankine-Hugoniot 


/(ml) = fiuR.) = /oo, 

and the Lax entropy condition 

/'(■Ul) > 0 > f'{uR). 

The above equation may be written in the equivalent form 

Ul > u> Ur, 


(2.3) 


(2.4) 


(2.5) 


where u is the unique state such that f{u) = 0. Integrating (2.1 1 ) in space over D, 
one obtains dt Jq udx = /(ur) — /(ur) = 0, which induces 


/ udx = / u^dx. 
In Jn 


( 2 . 6 ) 
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2.2. Discontinuous Galerkin formulation. The DG method consists in defin¬ 
ing a discrete weak formulation of problem (2.1). The domain is discretized with a 
uniform grid = UjgzKj with cells kj = Xj^i/2 = {j + \)h and 

/i > 0 the space step (see Figure 2.11. 

2.2.1. Numerical solution and Legendre polynomials. We look for approx¬ 
imate solutions in the function space of discontinuous polynomials 


— {vh € L^i^h) ■ € fl/i}, 


(2.7) 


where Vp{Kj) denotes the space of polynomials of degree at most p in the element Kj. 
The approximate solution to problem (2.1) is sought under the form 


Uh{x,t) = '^(l)j{x)Uj{t), \/x e Kj, Kj 0, 


( 2 . 8 ) 


1=0 


where Uj are the degrees of freedom (DOFs) in the element Kj. The subset {(jPj ,..., (j)j) 
constitutes a basis of restricted onto a given element. In this work we will use the 
Legendre polynomials Lo<fe<p- The basis functions in a given element Kj thus write 
(j>j{x) = Lk{2{x — Xj)/h) where Xj = {xjj. 1/2 + Xj_i/2)/2 denotes the center of the 
element. Orthogonality of the Legendre polynomials induces 


(j)^j{x)4>'j{x)dx = 2^ Vkj G ^h, 0 < k,l <p, 


(2.9) 


where Sk,i denotes the Kronecker symbol, and from (2.8) we obtain the following 
expression for the mean value of the numerical solution 


{Uh)j{t) 


1 

h 



Uh{x,t)dx = Uj{t), 


Vkj s ri/i, t > 0. 


( 2 . 10 ) 


Likewise, the properties Lfe(±l) = (±1)^ induce the following expressions for the 
left and right traces of the numerical solution at interfaces Xj- 1 - 1/2 of a given element: 


Uj+i/2(*) ■= w?i^+i/2.0 = 

1^0 

p 

Uj_i/2{t) ■= Uh{x^-l/2^t) Vt>0, 

1=0 


(2.11a) 

(2.11b) 


where denote the right and left traces of the quantity w at a given position (see 
Figure [2T| ). Finally, the derivatives of the Legendre polynomials may be defined from 
the recurrence relation [5] 


dsLk = dsLk -2 + {2k - l)Lk-i{s), k > 1, Lo{s) = l, Li{s) = s. (2.12) 

Multiplying the above relation with Li{s), integrating over [—1,1] and using or¬ 
thogonality, one obtains the recursion 


Nk,i = Nk-2,i + 2Sk-ij, Nqj = 0 , Nij 


26oj, 1 < k < p, 0 < I < p, ( 2 . 13 ) 
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h 



X 


Fig. 2.1. Mesh with definition of left and right traces at interfaces Xj^-j^/ 2 - 


and hence N/^ i = 1 — (—1)^+* if fc > I > 0 and N/^ i = 0 if 0 < fc < Z for the entries of 
the matrix N defined by 


Nkj = j Li{s)dsLkds, 0<k,l<p. 


(2.14) 


Note that we know the explicit projection of the derivatives of the Legendre 
polynomials into the Legendre basis. Indeed, applying the recurrence relation (2.121 
successively one obtains 


dsL2k = ^(4Z — 1)L2;-i(s), dsL2fc+i = ^(4Z + l)L2/(s), 

1=1 1=0 

which may be rewritten under the general form 

fc-i 

dsLk = MkgLi{s), 


(2.15) 


1=0 


where the entries of the matrix M are defined by Mkg = (2/ + 1) ^ — for k > I > 0 

and Mk ; = 0 for 0 < fc < Z. 


in space of problem (|2.1|) reads: find Uh in Vf such that 


2.2.2. Space discretization. The semi-discrete form of the DG discretization 

(2.16) 


/ VhdtUhdx- TZj{uh,Vh) = 0, Vu/i G V^. 


Kj 


together with the initial condition 

/ Vhix)uh{x,0)dx = / Vh{x)uo{x)dx, VvhGVj^. 


The discretization of the elementwise explicit residuals in (2.16) reads 


TZjiuh,Vh)= I f{uh)dxVhdx- f h{Uf^,ul)vhdS, 


(2.17) 


(2.18) 


where Zi : 17^ x > K denotes a numerical flux consistent with the physical flux: 
h{u,u) = f{u). We will consider the Godunov flux for the theoretical analysis of 
sections |4] and 15.21 




min{/(u) : v G K .ujj"]}, if < 

max{/(u) : v &[ul,ul]}, if >u^. 




(2.19) 
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^ (U,,uj 


u 


h 


fK)=fK) 



h 

f(u) 

fK) 

(u,u) 


fK) 

f(uj+f(u;)-f(0) 


(b) 


Fig. 2.2. Values of (a) the Godunov flux ji2. 1 9^ and (b) the Engquist-Osher flux Jl2.21]i) in the 
set of states for a convex scalar physical flux. 


The Godunov flux results from the solution of the Riemann problem for (2.1) with 
piecewise constant initial data consisting in states and at the left and right of 


the interface. The values of (2.19) in the set of states are depicted in Figure [2^ a). 

Comparisons will also be given in the numerical experiments for the local Lax- 
Friedrichs (LLF) flux: 


h(f. 




) = 




+ rK -O- 


( 2 . 20 ) 


min(u^ i^/T) max(u^ Engquist-Osher flux (see Figure 


2.2 


where the constant a is a stabilization parameter defined by a > maxjl f'(ri)| : 

b)): 

( 2 . 21 ) 




\f'iv)\dv, 


We recall that fluxes (2.19) to (2.21) are Lipschitz continuous and monotone 

du-h{a,b)>0, dy_+h{a,b) < 0, ya,bGfla, (2.22) 

where and d^+h denote the partial derivatives of h with respect to its first and 

second arguments. Numerical experiments tend to show that the effect of the flux 
on the quality of the approximation decreases as the polynomial degree p increases 
The analysis in section and numerical experiments in section [^will support 
this trend and confirm the relevance of the theoretical analysis with the Godunov flux 
in sections Handlll 

2.3. Steady-state scheme. Looking for steady-state solutions 


Uh{x) = ^ (j)jix)Up Vx e Kj, Kj € 


(2.23) 


1=0 


of the numerical scheme (2.16) and projecting the discrete scheme into the function 
basis, one obtains for all Kj in flh and 0 < k < p: 

n’;iuh) -TZjiuh^cj)^) 

= j f{uh)d^(t>''jdx-h{u:^^,u^^0+{-lfh{u~_^,u^_^). (2.24) 
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The local conservation property 


h{ 


u. 


i+5 i+ 




(2.25) 


is obtained from TZj{uh, 1^,) = 0 where denotes the indicator function of Kj. We 
are interested in discrete solutions of the numerical scheme such that 


lim Uh{x) = ul, lim Uh{x) = ur. (2.26) 

X—^—OO X—^OO 


Therefore, we obtain lim 2 ;_>±oo h{uf^ ,u^) = foo and the elementwise residuals in 
cells Kj in may be rewritten under the more convenient forms 




(2.27a) 


= 1 ,U+ i) - /oo = 0, 

J 2 J 2 

(2.27b) 


= [ f{uh)da,(l))dx - (1 - {-if) foo =0, 1 < fc < p. 

Jkj 

(2.27c) 


3. Steady discrete shock solutions for general convex fluxes. In this 
section, we present some preliminary results obtained for a convex flux / that will 
be used in the next sections. We assume that the exact shock position Xc is known 
and we set jc = 0 the cell index containing the shock without loss of generality. Cells 
j < 0 and j > 0 will be referred as to supersonic and subsonic cells, respectively. Our 
objective is here to determine the structure of the numerical solution in cells upstream 
and downstream of the shock position and to provide information about the solution 
in the shock cell. 


The analysis will mainly focus on the use of the Godunov flux (2.19). For a steady- 
state solution to hold, the conservation property = f^o, with a coercive and 

convex flux, has one of the following two solutions 


Uf^ = Ul and u'l > ur, 
^ Ul and ujj) = ur. 


(3.1a) 

(3.1b) 


Indeed, equation f(u) = foo bas two roots u = ul or u = ur (see section 2.1). 


Hence, to solve b(it^ ,u'f) = foo, it is suficient to consider the solutions of Riemann 
problems with initial conditions either uf such that f{uf) = foo and u'^ in V,a, or uf 


ifl, 


4 

or M. = Ur. 


and u'f such that /(«)() = foo- Consider the first initial condition, then Uf^ = ul 

+ 


Since h is continuous, we now consider successive situations when u^^ 
varies. Consider first uf = ul- The following situations (i) ul < > (b) ur < uf < 


Ul, and (iii) u'f < ur lead to the following expressions of the Godunov flux (2.19): 


(i) h{uL,uf) = f{uL), (ii) h(uL,ul) = f{uL), and (iii) h(uL,uf) = f(uf) > foo- 
The other situation, uf = ur, leads to h{uR,u'^) = mm{f{u),f{u^)} < foo when 


Ur < M + 


, where u is defined by (2.5), and h{uR,u'^) = f{u'^) > foo when ujj) < ur. 


We therefore obtain (3.1 1 ) and a similar step applies to infer ( |3.1| 3). 

In the following, we shall assume a stronger condition on the traces of the numer¬ 
ical solution at left and right faces of the shock cell: 


and > ur, 


'^\/2 ^ 


Ul and = '^R- 


(3.2a) 

(3.2b) 
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Indeed, the Godunov flux is no longer diff erent iable at points ) with 


> u such that /{Uf^) = /(u^) (see Figure 2.2'a)) and the strict inequalities 


in (3.2) ensure its differentiability in our analysis. We will give some comments on 


the implications of assumption (3.2) in section 3.3 We stress that the present results 


will show that the above assumption is not a restriction on the validity of our analysis. 


As a consequence of (3.2), the Godunov flux (2.19) reduces to a fully upwind flux 
at interfaces x±i/ 2 '- 




1 / 2 ) 




1 / 2 / 


(3.3a) 

(3.3b) 


From (3.3 1 ), the residuals in cell K-i no longer depend on DOFs in cell kq. Hence, 


the DOFs in cells / < 0 now become independent of DOFs in cells j > 0. Likewise, 


(3.3 d) implies that DOFs in cells j > 0 are independent of DOFs in cells j < 0. This 


is a key result for the theoretical analysis below and follows from the choice of the 
Godunov flux. 


3.1. Numerical solution in the subsonic and supersonic regions. The 

first result concerns existence of discrete shock profiles in cells j < 0 in the supersonic 
region when using the Godunov numerical flux. This result is local in the sense that 
it is based on the inverse function theorem in each cell and is valid in a neighborhood 


of the exact solution (2.2). The same analysis is then applied in the subsonic region. 


In the case of more general monotone numerical fluxes, we also give conditions that 
allow an exponential decay of perturbations away from the shock position. 

Theorem 3.1 (Solution in the supersonic region). Assume that the shock is 
strictly contained in cell jc = 0. Let p > 0, h > 0 and assume that the discrete 
solution satisfies Then in each cell j < 0, the constant discrete function in 

'Pp{Kj) define d by Wh ■ kj 5 x Wh{x) = ul is the unique solution to the numerical 
scheme (2.21) over a sufficiently small neighborhood Uj C 


Proof. We first note that, by consistency, Wh trivially satisfies the numerical 


scheme (|2.27|): TZj{wh) = 0, for all 0 < fc < p and j < 0. In other words, Wh is 


a stationary solution in the supersonic region. Let us prove that this is the unique 
solution in a sufficiently small neighborhood. Since ul > u (see Figure [ 2 )^ ), for each 
interface j + \ < — 5 , there exist neighborhoods Vj x Vj+i C x of the 

solution Wh over k, U where the numerical flux in (2.27^) reduces to the upwind 
flux: 


n%uh) = f{u-^,)- f^ = o. 


(3.4) 


U, = [U) 


j^0<l<p 


By assumption (3.2 j,), (3.4) holds also at interface X- 1 / 2 . Now let j < 0, define 
and let Gjj^i /2 '. —>■ be the application whose components 


are the elementwise residuals in (3.4) and (2.27 


= 0 ’ 0<k<p. 


(3.5) 


Setting = (ml, 0, ... ,0)^ in we have Gj+i/ 2 (UL) = 0 by consistency 

and it is sufficient to prove that the mapping Gj_|_i /2 is invertible over a neighborhood 
Wj C of Ul. The partial derivatives of Gj+ 1/2 take values 
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5c7iG°+i/2 = /'(Mj + 1 / 2 ) = f{uL), 0<1<P, 

f f {Uh)(t>’^jdx(t>)dx = f'{uL)Nk,h 0 <k<p, 0<l<p, 

J Kj 


^U]G^j + l/2 


at the point Ul, where the coefhcients Nk^i are defined by (2.14|. The application 
Gj+il 2 is continuously differentiable over and its Jacobian |Vu 3 Gj_|_i/ 2 l has the 
following form at the point U^: 


/'(ul) 

/'(ml) 

/'(ml) 

f'[uL) 

2/'(^l) 

0 

0 

0 

0 

2/'(ml) 

0 

0 


0 

2/'(ml) 

0 

il-{-lY)f{uL) 

{i + {-iY)riuL) 

il-{-lY)f{nL) 

0 


and reduces to iVu^Glu^ = {—2Y{f'(uL)Y^^ Y 0 because /'{ul) > 0. Applying the 
inverse function theorem, there exists a neighborhood Wj C of where the 

mapping Gj_|_i /2 is invertible. We conclude with Uj = Vj nw,n 

A similar result holds for the subsonic region that we formulate below and whose 
proof uses the same method as for Theorem |3.1[ 

Theorem 3.2 (Solution in the subsonic region). Assume that the shock is strictly 
contained in cell jc = 0. Let p > 0, h > 0 and assume that the discrete solution 
satisfies (3.2b). Then in each cell j > 0, the constant discrete function in Vp{Kj) 
defined by wt ■ kj 5 x Wh(x) = ur is the unique solution to the numerical scheme 


(2.21) over a sufficiently small neighborhood Uj C 


In the case of a general numerical flux, the upwind property (3.3) is no longer 


valid with the consequence that the DOFs in supersonic and subsonic regions are 
no longer uncoupled from other DOFs and oscillations may appear in cells close to 
the shock position. However, we now show that, under some assumptions on the 
numerical flux, these oscillations decay exponentially fast from the shock. This trend 
will be illustrated in the numerical experiments of section with the LLF flux and 
supports the relevance of the choice of the Godunov flux for the the oretic al analysis. 

Theorem 3.3. Suppose h : Ha x Dq — )■ K is a monotone (2.22), Lipschitz 


continuous numerical flux consistent with the physical flux f{u) and satisfies 


P ^ du+h{uL,UL) ^ ^ 
du-h{uL,UL) 


Q ^ du-hiuR,UR) ^ ^ 
du+h{uR,UR) 


(3.6) 


Assume that there exist J <0 and J'^ > 0 such that the discrete solution satisfies 
Uh{x) =ux + Y7i=o^]4>]{x) in cells j S {J“, J+} with l‘'jl - e < 1 and 

ux = Ur if j <0 or ux = Ur if j > 0. Then, the amplitude of oscillations of the 
numerical solution around the exact solution decays exponentially fast as |j| —>■ 00. 

Proof. We first consider supersonic cells j < J~ and proceed by induction. Let 
an interface j + 5 < J~ + ^ and assume that the property holds true in cell j + 1: 
Uh{x) =ul + ‘?j+i^j+i(a^) for all x S Kj+i with Yd^o l^j+il ^ define 

Gj+ 1/2 : X —)• RP+i the application whose components are the elementwise 


residuals in (2.27 j,), or equivalently (3.4), and (2.27:): 


/^k 

‘^i+ 1/2 


(U 2 , U,+i) := 'R'){uh) = 0, 0<k<p. 
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Note that Gj+ 1/2 also vanishes at (11^,11^) by consistency. By assumption of 
monotonicity of h over Oq x Oq, there exists a neighborhood Vj_|_i /2 C x 
of (UljUl) where Gj+ 1/2 is continuously differentiable and has partial derivatives 

+ + = i-^y9u+Hu~+i/2,U'^+l/2) = i-^y^u+HuL,UL), 0 < I < P, 


du‘ G^+1/2 


f f {Uh)4>\dx(t>’}dx = f{uL)Nk,i, 0 <k<p, 0<l <p, 

J Ka 


^u]^^G]+i/2 = 0. 0 < fc < p, 0 < Z < p, 


at (UljUl), where the coefficients are defined by (2.14|. In the same way as in 
proof of Theorem |3.H the Jacobian becomes 


|Vu,G, + i/2|(Ul,Ul) = {-2f'{uL)ydu-h[uL,UL) y 0 


from /'{ul) > 0, monotonicity of h and the first inequality in (3.6). Applying the 
implicit function theorem, there exists a neighborhood Wj_|_i /2 C about Ul and 

a continuously differentiable application p : Wj+ 1/2 3 U^+i i-A = (p(Uj+i) S 
such that Gj+i/ 2 (p(Uj+i),Uj+i) = 0. Setting Uh\^- = ul + ^ Taylor 

development of about U l reads 


^i+1/2 (U„ U,+i) = G°+i/ 2 (ul, ul) + h{uL, ul) 

1=0 

+du+h{uL,UL)'^{-iyyj +0{€‘^ul) 

1=0 

~ du-h{uL,UL){u~_^_^^^ - Ul) + i 9 „+/i(ul,ul)(m +^^/2 “ '^l)- 


We thus obtain that 


^j+ 1/2 


-UL^- 


du+h{uL,UL) ^ 

du-h{uL,UL) 


^j + 1/2 


- Ul). 


(3.7) 


From (3.6), the left and right traces of the numerical solution have the same sign. 
Likewise, a Taylor development of the applications G^_|_j^y 2 with 0 < k < p about Ul 


reads 


G^"+i/ 2 (U„ U,+i) = 0 = -f(uL) + Oi^ul), 0<k<p, 

1=0 

which induces <;j- = 0{e^u\) for 0 < / < p. The last coefficient is defined through the 
relation (3.7) and reads 


p ^ di^UL^ Ul) , _|_ 


S - - 


du-h{uL,UL 


■{u: 


'i+1/2 


- Ul). 


From (3.6), we have |i9„+/i(ul, /i(ul, ul)| < 1 and the last coefficient of 

the perturbation, y, thus decay exponentially fast as j —>■ —00, while other coefficients 
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are always smaller than . Therefore, the induction holds true in cell Kj and the proof 
is complete since it is assumed to be true in cell J~. 

The proof is similar for subsonic cells. Introducing the application 


:= =0, 0 < fc < p, 


defined by (2.27 d) and (2.27:) for j > J+ and assuming small perturbations about 
, 0)^ in the implicit function theorem and a Taylor development 

r|j) for 0 < A: < p and 


about (UfljUfl) give = 0{e^ 


p du-h{uR,UR) _ 

S - -IT—^7"^—(u,-1/2 - Ur). 

du+HUR,UR) 


□ 


Note that the assumption (3.6) holds for the LLF flux. Indeed, we have du-h{u , 
M+) = {f'{u~) + q)I2 and du+h{u~,u^) = — q)/2 with —a < f'{uR) < 

0 < /'{ur) < a. Numerical experiments of section 6.1 will illustrate the results of 
Theorem 3.3 In [32], Smyrlis used similar arguments to demonstrate the exponential 
decay of oscillations in the vicinity of a stationary shock for the Lax-Wendroff scheme 
observed in precedent numerical experiments m- In this case, monotonicity is lost 
and oscillations change sign from one cell to another. 

The proof of Theorem |3.3| shows that the oscillations are transmitted from a cell 
to its neighboring cell Kj by the numerical fluxes and that only the highest DOF is 
affected. In particular, setting UJ = 0 one recovers the uniform solution up to first 
order: = ux + This may motivate limiting techniques of the solution 

to suppress spurious oscillations and a first attempt in this direction will be shown in 
section 16.31 


3.2. Numerical solution in the shock region with the Godunov flux. In 

the case of the Godunov flux, the mean value of the numerical solution in the shock 
cell is explicitly known. This is the object of the following lemma. 

Lemma 3.4 (Global conservation). Let h> 0 and p > 0 and use the Godunov 


flux {2.19), then under the assumptions {3.2^^ the first DOF in cell kq is defined by 

(1 + Sc)ul + (1 - Sc)uR 


[/“ = {u)o = 


(3.8) 


where (u)o denotes the mean value (2.10) of the exact solution {2.2) in cell kq, and 


Sc = J^{Xc - xo) 


(3.9) 


is the relative shock position in kq and has values Sc = ±1 when the shock is at 
interfaces x^^ji. 

Proof. Setting vt = Ik^ into (2.16) and summing up the results over all cells kj 
in fl/i, one obtains 


4( S - 


t/0)= ^ 7^° = - ^ h^^^-h^_^=f{uL)-f{uR)=0, 


according to (2.26). Integrating (2.1 1 ) in space over D, one obtains dt f^udx = 
/(wl) — /(ur) = 0. Subtracting this result to the above equation, one obtains 


d 

dt 


('• i: 


KjG^h 


I u 

Jn 


{x)dx^ = h^ (ul^ - (■u)o) = 0 , 
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from Theorems 3.1 and 3.2 which gives Uq = {u)o from (2.17). □ 


3.3. Some comments on assumption 


The above results with the Go¬ 
dunov flux are based on assumption (|3.2|) on the traces of the numerical solution at 


interfaces of the shock cell. However, from (3.1) we observe that this assumption 
may be violated when either (i) < ul and or (ii) = ul and 

^ '^R- Our objective is here to show that, up to a shift in the index of the shock 


cell, jc, (3.3) still holds. This is a consequence of the following result: for oscillations 


in cell koj situations (i) and (ii) reduce to 


^- 3/2 “ '“- 1/2 


= UL, 


^- 1/2 


^ 1/2 RRi 


^- 1/2 


^1/2 ^1/2 ^3/2 U-iJ, 


1/2 ~ “ 3/2 


(3.10a) 

(3.10b) 


respectively. To prove these results, consider the local residuals (2.18) at steady-state. 


Setting Vh = Uh und imposing the conservation property (2.25), one obtains 

TljiUh, Uh) = g{u~+i/2) - 5(u+_i/ 2) - /oo(uJ+i/2 - ^j^- 1 / 2 ) 

= (/(O) - /oo) (M7+1/2 - u^-1/2) = ^ 


(3.11) 


where g{u) = f{v)dv and min(u^_j^^ 2 J ^^+ 1 / 2 ) — from the 

mean value theorem. The above equation has two solutions: either u7^j^y2 “ '“j"_i/ 2 ’ 
or /(Ci) = foe and min(u+_^/ 2 > "“ 7 + 1 / 2 ) < 0 

Now, suppose for instance that case (i) holds so that u~f ^^2 = '“R- Applying 


(3.11) with j = 0, we obtain either the trivial solution = ur, or /(^o) = foe 
with ^0 Ur- In the latter case, we have ur < ^0 = ul < 117/2 coercivity and 
strict convexity of / (see section 2.1). This latter result ca nnot hold because of the 


definition of the Godunov flux which imposes < ur by (3.1), so = ur. 

Then, every solutions such that uZi /2 ^ correspond to situations where 
the oscillations are in cell Kj^ with < —1. Indeed, by (2.19) one would have 


h('n_i/ 2 ’^^ 1 / 2 ) ~ f{uf^^ 2 ) und since h(itj^^ 2 > 1 I 1 / 2 ) “ /('“i/ 2 )> there exists a neigh¬ 
borhood of the uniform solution Uh = ur in cells j > 0 where it is the unique solution 
from The orem Hence, by setting kq the cell containing osci llations accord ing t o 


Theorems 


3.1 


and 


3.2 


only ii _]^/2 “ violates assumption (3.2). Using again (3.11) 
iZi /2 — ul imposes — ul- Similar arguments hold for proving that (3.10 2 ) 


corresponds to the only situation of case (ii) violating (3.2). 

To sum up, it is convenient to introduce an index of faces of the shock cell 


ic + k = Then, the situations violating (3.2) reduce to one of the following cases 


-2 

“t-l /2 “ '“ic + l /2 


= Ur, 


■“^ + 1/2 + 


= UR, 


-F I =±I 

^2 2 ■ 


Therefore, the Godunov flux is not differentiable at Xi^+if 2 only, /ii^+i /2 = 
j^if 2 )'- solution satisfies (|3.3|) and the analysis of section ^ remains valid, 
but the other one corresponds to oscillations in a neighboring cell. In the following. 


we will only consider solutions for which (3.3) holds without loss of generality. Among 


these solutions, the situations violating (3.2) will be analyzed in the linear stability 


analysis of section and numerical experiments of section 
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4. Steady discrete shock solutions for the Burgers equation. We are now 

interested in steady-state solutions of the invis cid Burgers equation, i.e., (2.1) with 
f{u) = \v?. According to Theorems 


3.1 


and 


3.2 


the uniform numerical solutions 


Uh = ul in the supersonic region and Uh = un in the subsonic region hold with the 
Godunov flux. We thus restrict our analysis to the solution of the numerical scheme 
in the cell containing the shock. The result now depends on the order of the numerical 
scheme and we give the solutions for 1 < p < 3 in the following theorem. 

Theorem 4.1. Consider the discrete DG scheme with the Godunov flux (2.19). 
Let the exact position of the shock he in cell Kq and assume that the internal traces 
satisfy (3.2). Then, for h > 0, according to the polynomial degree, p < 3, and the 
relative position of the shock in the cell, Sc, defined by (3.9), the solution of the discrete 
scheme (2.21) reads: 


p 

Uhix) = ul'^ui(I)q{x), \/x G ko, 
1=0 


with uq = Sc and if p = 1: 


ifp = 2: 


- 1 < Sc < 1, Ul = -\/3(l - s2); 


(4.1) 


ifp = 3; 


Sc G Di, 


Sc G I?2i 


1 < Sc < 

-1 \ 

e s 

II II 

0, 

\/5(l - s2), 

(4.2a) 

2 

-F < Sc 

4 \ 

1 Ul = 


(4.2b) 

O 


[ U2 = 

-|sc, 


2 

3<sc 

<1, ] 

f Ul = 

1 U2 = 

0, 

-\/5(l - s2); 

(4.2c) 


1 /c 

5 V 2 V ^ 

- 54s2, 



U2 = 

-7sc, 



(4.3a) 

U3 = 

-Ul, 




Ul = 

-1 r 

20\/ 3sc A 2 

- 7(17s2 - 

- 6)3-L 49sc(17s2 - 6)2A3^^ 

-L 


26s2(17s2 - 

OAf- 

Sc(463s2 - 246)A3- 



26ScA2Ay^ -1- ScAg^^ 

^ + a“iA 


U2 = 

12 [ 

7(17s^-6) 
^1/3 
_ 3 

+Ar], 

(4.3b) 

U3 = 




4207sc(17s2 

:-6)2A3/'- 420s2(17s2-6)A3/% 



A2(2sc(5969s 2 - 2889) -L 420s2Ay% 
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Sc € 2 ? 3 , 


Ui = 


,0^ (108252 - 386051s2 + 3356194) + 


U2 

U3 


V3Ai (ll64 - 24035^ + aI^^{26sc 
a 44588 - 9144 + 10895^) + 
A3^4703 - 11254) 

J_r_7„ _ 7{17sI-6) 

12 L 


7A 


1/3 


-, 1/2 


/W 


-A 


1/31 


(4.3c) 


7ui 


-^[35(12-1614 + 2364)+ 


108 ( 1 -s 2 )A(' 

STSA^Ag/^ - 7) + 5A3 /^s(79 - 1004) + 

2A^/4324-17)], 


^2 = -D: 


= (-l>-\/^)U(i,l); 


A 3 = 7 A 2 , 

(4.4a) 

(4.4b) 


where Vi = i), 

A 3 = — 7 A 2 , A 2 = —A 2 and 

A2 = 343Ai - 419s3 + 279sc, 

Ai = [7776s® - 10008s^ + 33595 ^ - 56] 

Proof. See Appendix]^ □ 

Several conclusions may be inferred from Theorem |4.1[ First, for a polynomial 
degree p, ie., a given DG scheme, the solutions are parametrized by the relative 
shock position in the cell only. In particular, the numerical solution is independent 
of the mesh size and the amplitude of the oscillations are proportional to the shock 
strength ul — un. These oscillations occur for p > 1 and vanish if and only if the 
shock is at an interface: Sc = ±1. Then, the coefhcients are displayed in Figure 4.1 a.- 
c) and present continuous evolution with Sc- The validity of the assumption 
is illustra ted in Figure j+^d) and is satisfied according to the ranges of solutions in 
Theorem |4.ll We observe that for p = 2 (resp. p = 3), the positions Sc = ±i (resp. 


Sc = ±g) correspond to situations where = '“ 1/2 = or ur (see section 3.3). 

It is remarkable that assumption (3.2) selects the solutions over the whole range Sc in 


[—1,1]. Finally, the oscillatory behavior of the solution may lead to entropy violating 
solutions where the sign of eigenvalues changes locally upstream and/or downstream 
of the shock position. Since the numerical solution can potentially change sign p times 
in a cell, this may lead to nonphysical situations. Numerical experiments in section 
l^will show that those situations may exist for p > 2. This feature is different from 
what has been shown in [18] for the high-order RBC schemes where the oscillations 
were seen to monotonically decrease for a shock position moving from the center of 
the cell to its faces. 

5. Linear stability of steady shock profiles. 

5.1. Linearized operator. In this section, we are interested in the linear sta¬ 


bility of the DG scheme around steady shock solutions Uh satisfying (2.27). We will 


mainly focus on sufficient conditions for instability of the numerical scheme. These 
conditions will then be used to analyze non-convergence to steady-state of the DG 
scheme observed under certain conditions in the experiments of section 

To this end, we first consider a forward Euler method for the time discretization, 


setting w/j = in (2.16), we obtain 


U. 


k{r 


•+1) _ jjk(r, 


+ Afc7^f(u40, VjeZ, 0<k<p, 


(5.1) 
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Fig. 4.1. From (a) t o (c) : 
cell K,o given by Theorem 4-i 


evolution of coefficients in the function basis of the solution in shock 
In (d): evolution of the quantities + ul (open symbols) and 


^ 1/2 symbols) with = —nji = 1 


where = Uj{nAt) is a function of time, Xk = (2fc + 1)A and A = ^ with At > 0 
the time step. We shall assume that the numerical flux is differentiable at every point 
Xj+i /2 with respect to the left and right traces steady solution. 


Let (j) : V] 


■UP. y('l) 


!-)■ be the map defined by (5.1) where 

is defined by (2.23) with components Now define the Gateaux derivative of (p 

at point Ufi in the direction ui/i by L : —>• Vj^;wh >—t Lwh = d4>{uh]U)h)- Using the 

vector notation Lw for the components in of Lwh, the Gateaux derivative of the 
DG operator (5.1) along the direction Wh reads 


(Lw)^^ 


Wf + Xk lim 

■' e-fO 


n'^iuh + ewh) -n^jiuh) 


0 < k <p. 


e 
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Using expression (2.24) for the local residuals, we obtain 


TZ’^{uh + ewh) = j f{uh + ewh)dj;(l)jdx - 


+(—l)^h(M 1 + ( 


ti)- 


Substracting TZ^{uh), dividing by e, and letting e tend to zero, the components 
of Lw read 

(Lw)J = Wj +\k[^ f Whf{uh)d^(l)^dx - 
J Kj 

+ {-^)''iw~_idu-hj_i +wt i5„+hj_i)). (5.2) 


The linearized operator L is thus made of three diagonals of blocks of size (p 
1) X (p + 1) with entries 




(5.3a) 


(Lj,i)fc,i = + f {uh)4''‘jdx4>'jdx - - (-l)'=9,,+ /i^_ i),(5.3b) 

Jkj 

^kdu+hj_^_i^ (5.3c) 

for 0 < k,l < p and j in Z. 

For the shock profile Uh to be linearly stable, it is necessary that the spectrum of 
L contains only eigenvalues p with modulus lower than unity, |p| < 1, and semisimple 
eigenvalues with unit modulus, |p| = 1 |3]. In the following, we will focus on sufficient 
conditions for instability of shock profiles. 

5.2. The case of the Godunov numerical flux. The upwind character of 
the Godunov flux allows to specify the eigenvalues of the linearized operator in the 
following proposition. 

Proposition 5.1. Under the assumptions of Lemma \3.4\ the spectrum of the 


linearized operator (5.3) for the Godunov flux (2.19) reduces to the spectra of the 


following matrices of size p + 1; 

(L‘j<o)k,i = Sk,i + Xkf {uL){Nk^i — 1 ), 

(Lo)a:,/ — T Xk j f {uh)fodx(fQdx^ 

J Ko 

m>o)k,l = 4.1 + Afe/'(ufl)(iVfc,; + (-!)'=+'), 


(5.4a) 

(5.4b) 

(5.4c) 


for 0 < k,l < p, where Uh denotes the solution of 1(2.2'/^ and Nkj is defined by ( 2.14 )■ 
Proof. Under assumption (3.2), the Godunov flux reduces to the upwind flux 


(3.3). Therefore, for j 0 the integral in (5.3 d) reads 


/ f {uh)(t>\dx(t)^^dx = f{ux)Nu^u 

with ux = ul for j < 0 and ux = ur for j > 0. After simplification and omitting 
the double subscript for diagonal blocks, the linearized operator reduces to 


Xj,]-i)k,i = (—1)^A/c/'('Ul), (Lj)fc,/ = 4,; + Xkf'{uL){Nk,i — 1), = 0, j < 0, 
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T able 5.1 

Eigenvalues of the iteration matrix j'5. ,?[ ) associated to the DG scheme with the Godunov nu¬ 
merical flux and a forward Euler time integration for the discretization of the Burgers equation. 

- 2 1 7 5 

Here, A/. = \ f'{uif), Xu = Xf'{uji), X = uj^X, 71 = 33 —33, and 72 = 36 + 36 . 


p 


j < 0 

J = 0 

j >0 

0 

(-1,1) 

1 — Al 

1 

1 ± A r 

1 

(-1,1) 

1 - Ai(2±f7/2) 

1, 1-2AV3(1-s2) 

1 + Xr{2±iV2) 


(-1,-i) 

1 — Al( 3 -f 7i), 

1, l±2fA X ... 

1 + A_r( 3 ± 7i), 


^-3scv'5(l - s2) - 6(1 - si) 

2 

(-14) 

1-^(6-7i 

1, 1-AV3(4-9s2), 

1 - 2Aa/ 3(4 - 9s2) 

1 + 4^(6 -71 


(1,1) 

±72*) 

1, l±2fA X ... 

± 72 *) 


^3sc\/5(l - si) - 6(1 - si) 


= (-l)"Afc/'K), (Lo)fc7 = 4,z + Xk f nuh)4d,4dx, 


(Lo,i)fe,i = and 

= 0, (L,)fc7 = 5k,i+\kf\uR){N^g + (-1)'=+'), = -{-iy\kf{uR), 

for j > 0. Therefore, L is lower block triangular for j < 0 and upper block triangular 
for j > 0. The eigenvalues of the matrix are thus the eigenvalues of the diagonal 


blocks which are constant and equal to (5.4). □ 


5.3. Application to the Burgers equation. As an application, Table[5TT|gives 


the eigenvalues of the three different blocks in (5.4) for polynomial approximations 


0 < p < 2 in the case of the Burgers equation f{u) = We stress that eigenvalues 
in blocks j 0 remain valid for a general physical flux /. 

Stability in the shock cell imposes A = urX < (3(1 — fo r p = 1, A < (3(4 — 
for p = 2 and |sc| < |, but the linearized operator (5.3) is unconditionally 


unstable for p = 2 and |sc| > §• The DG method is indeed unstable for p > 1 and a 
forward Euler method. These modes may be s tabilized by using Runge-Kutta schemes 
of sufficient order (see below and section 6.2). For p = 2 and Sc £ (—1, — |), the real 


part of eigenvectors associated to the unstable eigenvalues has r = ( 0 , 0 , 2-y5(l — 3 ^) + 
5sc)^ as components in the basis of 'P 2 {ko). For p = 2 and Sc € (|, 1), their real part 
reads r = (0,0, 2-y5(l — s^) — 5sc)^. They only affect the highest DOF and reach 
largest values at faces of the shock cell. 


Let us consider the linear stability of situations violating assumption (3.2) ac¬ 
cording to section |3.3[ For these points. Table |5.2| displays the eigenvalues and 
eigenvectors in the shock cell Lq. We observe that for Sc = if and 1 < p < 2 , 
p = 1 is a non semisimple eigenvalue of L. Indeed, A = ^ > 0 hence p = 1 is 
not an eigenvalue of blocks Lj<o or Lj>o (see Table 5.1), but it is an eigenvalue of 
Lq with algebraic multiplicity p -I- 1 and geometric multiplicity of 1 as indicated in 
Table [5(2} This property holds for p = 3 because eigenvalues of Lj<o and Lj>o satisfy 
p 7 ^ 1. Indeed, they are of the form 1 — A^p/ and 1 -I- XrHi, respectively, with p; £ 

{-4+(73-4)1/2±^+, _4-(y3-4)i/2±^-}, ^3 = 102/3(*y6-2)-i/3 + (i0(iv^-2))i/3 
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Table 5.2 

Eigenvalues and eigenvectors in the shock cell of the iteration matrix associated to the DG 
scheme with the Godunov numerical flux and a forward Euler time integration for the discretization 
of the Burgers equation. Here, X = u^X 


p 

Sc 

eigenvalues 

eigenvectors 

1 

±1 

{1,1} 


2 

±1 

{1,1,1} 

/O 0 0\ 
0 0 0 
Vi 0 oj 


=b^ 

3 

{1,1,1} 

(11 
\1 0 0/ 


3 

±1 

{1,1,1,!} 


/O 0 0 0 
0 0 0 0 
0 0 0 0 
1^1 0 0 0 

\ 

/ 


{1,1,1,1-f} 


1^^^ 0 0 0 \ 

^ T# 0 0 

0 1 0 

\ 1 0 0 1 j 



and 7 ^ = i{8 + 73 ± 8(73 — 4)“^/^)^/^. Note that the eigenvector associated to /r = 1 
has components along the highest DOF only. This property will be used in section [63 
for stabilizing the DG scheme. 

Table [53] also gives the eigenvalues and eigenvectors for shock positions — 1 < 


< 1 where the strict inequality in assumption (3.2) is violated. This situation 

Though different from 


occurs for p = 2 and p = 3 as highlighted in Theorem [4.1 
results for Sc = ± 1 , these results show that p = 1 is also a non semisimple eigenvalue 
and the DG scheme will be unstable. At this time, a general result about instability 
has not been achieved and is beyond the scope of the present study. 

Stability of the DG scheme requires high-order Runge-Kutta schemes [S]. How¬ 
ever, the scheme may remain linearly unstable at points violating assumption (3.2). 


As an example, consider the second-order and strong stability preserving Heun scheme 
whose linearized operator reads Lrk 2 = L -f i(L — I)^. Solutions of equation 
^J'RK 2 = M ~ 1)^ = 1 are p = ±1. We note that L is lower triangular by 

blocks for rows j < 0 and upper triangular by blocks for rows j > 0, so is 'Lrk 2 - As 
a consequence, the stability of 'Lrk2 reduces to the stability of its diagonal blocks. 
Moreover, the transformation from L to 'Lrk2 results on the same operation of its 
diagonal blocks. Now, we note that p = ±1 is not a root of plrk 2 = 1 for blocks Lj<o 
or Lj>o for p < 3. Indeed, according to Table [5.1 j and precedent remarks, the only 
roots are Xr = A/'(ml) = 2 for p = 0 and Xr = 2/(3 -I- 71 ) ~ 0.5498 for p = 2 which 
do not satisfy the usual GFL condition max{|AL|, |A_r|} < l/(2p + 1) [S]. The same 
analysis holds for Xr = Xf'{uR). 

6. Numerical experiments. We consider the inviscid Burgers equation, f{u) = 

over D = [0,1] with boundary conditions m( 0) = —w(l) = 1. The numerical so¬ 
lution is obtained by using a method of lines. The semi-discrete equation (2.16) is 


advanced in time by means of an explicit third-order and strong stability preserving 
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Runge-Kutta method [3T]. We look for steady-state solutions of (2.16) of the form 
Uh{x) = lim„_>,oo Uh{x,nAt). The time step is set at 


At = CFL X min 




\ frf M ■ ^0 ^ ! 


where CFL = l/(2p + 1) according to [5] and the maximum eigenvalue is evaluated 
at quadrature points of the element Kj. 

As suggested in [TB], the final shock position for the Burgers equation is set 
through the following initial condition 

, ^ f l — 2(l — u)x, if a; < i, 

Uo(a:) - I 2 u + l- 2(u + l)a:, 


with —2<u<2. Using (2.6), we obtain Xc = 1/2 -I- u/4. 


We note that the evaluation of the volume integral in (2.18) is done by using Gauss 


quadrature which may be inexact for a nonlinear flux. The present results have been 
obtained by using a numerical quadrature of sufhcient order to integrate it exactly in 
the case of the Burgers equation: p -|- 1 Gauss-Legendre points are used for p < 2 and 
p -|- 2 points are used for p = 3. The extra point for p = 3 allows to quantitatively 
compare the solution of the numerical scheme obtained from the theoretical analysis 
in Theorem 14.11 with the solution obtained from a numerical calculation. 


6.1. Structure of steady shock profiles. Figures [6T] to |6.3| display the steady- 
state solutions to the DG scheme obtained with the Godunov flux ( |2.19[ ) for nine 
different shock positions (3.9). We compare solutions obtained from the theoretical 


analysis in Theorem |4 .1 1 with the solution obtained from a numerical calculation with 


IV = 20 cells and jc = 11. Solutions are also compared to the exact solution (2.2). 


The solution remains uniform in the supersonic and subsonic regions x < x_i /2 
and X > Xi /2 whatever the polynomial degree which confirms the conclusions from 


Theorems 3.1 and 3.2 However, for p = 2 and Sc = 1 or for p = 3 and Sc = ±1, the 


solution from the numerical calculation appears to be oscillatory outside the shock cell 
and differs from the expected predictions. In all these situations, the calculations did 
not succeed in converging to a steady state and therefore do not satisfy the discrete 


scheme (2.27). The analysis in section 5.3 predicts instability in situations Sc = ±1 
for 1 < p < 3, Sc = ±1 for p = 2, and Sc = ±g for p = 3. The analysis was limited to 
a second-order Runge-Kutta scheme, but results of section [R2] show that it holds for 
a third-order scheme. We recall that assumption ( |3.2[ ) may be violated for p > 1 at 
these points (see Theorem |4. 1 1 and Appendix [A| where the God unov flux admits two 


equal values /i_i /2 = or hi /2 = fiui/ 2 ) (see Figure 4.1 'd)). Both values may 


be solutions of the numerical scheme as long as that the conservation of the scheme 


(2.25) is respected. Our numerical experiments tend to show that, when converging 


to the steady state, the Godunov flux changes for a solution to another following a 
cyclic pattern with the consequence that the flux balance in the cell is periodically 
modified (see Figure [U^ b)). In some cases, e.g., p = 1 and Sc = ±1 or p = 2 and 
Sc = — 1, the convergence to steady state was reached but at very low speed. 

The solution oscillates in the shock region. The theoretical solutions of The- 
4.1 agree very well with the numerical calculations. The oscillations present 


orem 


amplitude lower than two times the shock strength ul — ur. Likewise, for p = 2 and 
p = 3 the oscillations of the polynomial solution in the shock cell may lead to non¬ 
physical solutions where the sign of the eigenvalues f'{uh) = Uh changes sign locally 
more than one time. 
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Uj, (theory) 

1.5 

U(, (theory^ 

1.5 

(theory^ 



0.5 

0 

■0.5 



0.5 

0 

■0.5 




-1.5 
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Fig. 6.1. Steady-state solutions in the shock cell kq for p = 1 and the Godunov flux ^77^. 


In Figures 6.4 and 6.5 we present numerical experiments with the LLF flux (2.201 
and Engquist-Osher flux (2.21). The results in the shock cell are remarkably similar to 
the theoretical solution for the Godunov flux (2.19). This observation is in agreement 
with precedent numerical evidence of comparable resolution of the DG scheme with 
different numerical fluxes when the polynomial degree p is increased [8l[28]. Note that 
the Engquist-Osher flux reduces to the Godunov flux as soon as the left and right 
traces satisfy u~ > it and < u, where u is defined in section 2.1 For instance, it 
may be easily checked that this holds for |sc| < ^ for p = 1 and |sc| < ^ for p = 2 
from Lemma 3.4 and Theorem 4.1 The main difference consists in the neighboring 


cells of the shock region where oscillations occur with the LLF flux. These oscillations 
are a consequence of Theorem |3.3| which details the mechanism of transmission of 
oscillations at interfaces with a monotone numerical flux. As an illustration, Table 


6.1 gives the values of the left and right traces at interfaces Xj^i /2 of the oscillations of 
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Fig. 6.2. Steady-state solutions in the shock cell kq for p = 2 and the Godunov flux ^77^. 


the numerical solution for 8 < j < 10 in the supersonic region. The trace ~ 

for j = 10 co rresponds to the left trace in the shock cell jc = 11. As expected by 
Theorem 


3.3 


the sign of the jumps ~ ul ts conserved through the interfaces 

as a consequence of the monotonicity of the LLF flux. The last column provides the 
values of DOFs oi Uh — up in each cell scaled by the jump — up. The values 

are always quite lower than u^+ 1/2 ~ every DOF for 0 < Z < p is several order 

of magnitude lower than the last DOF I = p. 


6.2. Linear stability of steady shock profiles. Figure |677| presents the spec¬ 


tra of the linearized operator (5.3) around steady-state solutions for different polyno¬ 


mial approximations and numerical fluxes when varying the relative shock position in 
the range [0,1]. Results are given for explicit first- and third-order time integration 
schemes. We are here interested in the slow convergence and even non convergence 
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Fig. 6.3. Steady-state solutions in the shock cell kq for p = 3 and the Godunov flux ^77^. 


observed from the numerical experiments in the preceding section and illustrated in 
Figure [ 6 ^ Considering only p = 1 and p = 2 approximations thus appears sufficient 
for our purpose. The results for the Godunov flux have already been determined from 
the stability analysis in section and will be used for comparison. 

For the second-order approximation, p = 1, every eigenvalue is contained in the 
unit disc. In the case of the Godunov flux, the pairs of complex conjugate eigenvalues 
correspond to modes located in cells J 7 ^ 0 and associated to uniform flow, while 
real eigenvalues correspond to modes in the shock cell as indicated in Table |5.1[ 
The spectra obtained for other numerical fluxes look similar even if a scattering of 
eigenvalues in cells j 7 ^ 0 is observed for the LLF flux. Using a third-order time 
integration scheme is seen to lower the modulus of eigenvalues as expected. 

The third-order approximation, p = 2, with the forward-Euler scheme exhibits 
unstable modes when the shock position tends to an interface. Here again, the spectra 
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(a) sc = - ^ (b) Sc = 0 (c) Sc = I 

Fig. 6.4. Steady-state solutions in the shock cell kq obtained with the LLF flux (2.20]) forp = 1 
(top row), p = 2 (second row), and p = 3 (bottom row). 


for every numerical scheme look similar. The complex eigenvalues are concentrated 
in the right half-plane, while some real eigenvalues are negative and may become 
unstable. From the stability analysis for the Godunov flux in Table |5.1[ one may 
distinguish among the former modes between stable eigenvalues associated to the 
uniform flow and unstable eigenvalues associated to the shock. For other numerical 
fluxes, one recovers the eigenvalues in cells j ^ 0 concentrated at the same locations 
1 — Al(3 -I- 7i) ~ 0.272 and 1 — ^(6 — 71 ± 721) — 0.464 ± 0.610i. Likewise, modes 
associated to the shock are also very similar. The eigenvalues are concentrated at the 
unit circle and become unstable as the shock position tends to the interface. We stress 
that lowering the CFL value (see eigenspectra for p = 2 and Sc = 19/20) sho ws that 
these modes remain unstable as predicted in Table 5.1 for |sc| > §• Figure 6.8 displays 


the structure of associated eigenvectors for two different shock positions corresponding 
to unstable modes. In section |5.2[ it was shown that these modes only affect the 
highest DOF and result in a quadratic evolution with support in the shock cell. The 
structure of the modes for the Engquist-Osher and LLF fluxes is very similar and 
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Fig. 6.5. Steady-state solutions in the shock cell kq obtained with the Engquist-Osher flux 
\2.21^ for p = 1 (top row), p = 2 (second row), and p = 3 (bottom row). 
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Fig. 6.6. Convergence histories to steady-state solutions of the Burgers equation for different 
kind of numerical fluxes (GOD: Godunov, LLF, OSH: Engquist-Osher) and relative shock positions 
0 < Sc < 1 with p = 2 and N = 20. 
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Table 6.1 

Illustration of Theorem \3.3\ for the LLF flux. Results obtained in the supersonic region for 
1 < p < 3 and Sc = —| (top), Sc = 0 (middle), and Sc = | (bottom). 


p 

j 

^i+i/2 

U++1/2 - UL 


{Uj —ulSi^o)o<i<p 




10 

-1.141e-02 

-2.245e-01 

9.639e-05 5.075e-02 

1 

9 

3.215e-05 

1.137e-02 

1.515e-08 2.827e-03 


8 

-2.585e-10 

-3.215e-05 

3.452e-12 8.038e-06 


10 

7.513e-02 

6.335e-01 

9.049e-04 1.425e-13 1.194e-01 

2 

9 

1.359e-03 

7.513e-02 

2.461e-06 1.198e-12 1.809e-02 


8 

4.618e-07 

1.359e-03 

2.433e-ll 7.067e-ll 3.396e-04 


10 

1.632e-01 

1.003e-00 

1.918e-03 4.649e-05 l.lOOe-03 1.635e-01 

3 

9 

-6.305e-03 

-1.648e-01 

1.721e-05 1.595e-08 9.839e-06 3.823e-02 


8 

9.900e-06 

6.303e-03 

1.103e-09 9.581e-12 6.455e-10 1.570e-03 


10 

9.637e-02 

7.320e-01 

2.187e-03 1.338e-01 

1 

9 

-2.365e-03 

-9.957e-02 

9.357e-06 2.374e-02 


8 

1.395e-06 

2.363e-03 

1.362e-10 5.901e-04 


10 

9.637e-02 

7.320e-01 

1.294e-03 1.293e-13 1.329e-01 

2 

9 

2.214e-03 

9.637e-02 

5.089e-06 9.095e-13 2.297e-02 


8 

1.224e-06 

2.214e-03 

4.613e-ll 4.139e-ll 5.529e-04 


10 

-6.651e-02 

-5.736e-01 

5.488e-04 5.358e-06 3.137e-04 1.157e-01 

3 

9 

1.061e-03 

6.624e-02 

1.215e-06 1.888e-10 6.946e-07 1.602e-02 


8 

-2.817e-07 

-1.061e-03 

3.889e-ll 5.812e-ll 6.699e-ll 2.653e-04 


10 

1.581e-01 

9.842e-01 

4.483e-03 1.651e-01 

1 

9 

-6.460e-03 

-1.6695e-01 

4.157e-05 3.865e-02 


8 

1.035e-05 

6.446e-03 

2.773e-09 1.606e-03 


10 

-2.489e-01 

-1.182e-00 

5.007e-03 6.805e-14 2.054e-01 

2 

9 

-1.390e-02 

-2.489e-01 

7.744e-05 3.599e-13 5.578e-02 


8 

-4.799e-05 

-1.390e-02 

1.656e-08 6.728e-12 3.451e-03 


10 

-6.6171e-03 

-1.690e-01 

1.849e-05 1.797e-08 1.056e-05 3.913e-02 

3 

9 

1.0901e-05 

6.614e-03 

1.276e-09 8.861e-12 7.436e-10 1.648e-03 


8 

-2.9690e-ll 

-1.090e-05 

4.032e-09 5.432e-09 6.223e-09 2.739e-06 


mainly concentrated in the shock cell. These modes are damped with the Runge- 
Kutta scheme, but are clustered around the unit circle and remain unstable for large 
Sc values. These modes are expected to slow-down or prevent the convergence of 
some computations as observed in Figure [Rbj i and in the numerical experiments in 
the preceding section. 


6.3. Stabilization with the spectral vanishing viscosity method. The an¬ 
alytical results in Theorem |3.3| and section [5?3| and numerical experiments in Table [6T] 
suggest the possibility of removing oscillations in the supersonic and subsonic regions 
by damping the higest DOF, say t/J with j ^ 0. Likewise, the structure of solutions 
in Theorem 4.1 and the stability analysis in section (see, Table 5.21 support 
the necessity of acting on a larger range of DOFs, say {Uj)mj<i<p where mj > 0, 
whenever the solution becomes discontinuous in cell Kj. These observations motivate 
the application of the spectral viscosity method [22J [23j to the DG discretization. In 
this method, one supplements the explicit residuals (2.18) with an artificial viscosity 
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of the form 


—e 


Q dx 


where e > 0 plays the role of a viscosity coefficient and Q denotes the spectral viscosity 
operator 


Qvh = ^ yvh = Kj G Qh, 


1=0 


1=0 


with Qj = exp ( — if I > rUj and Qj = 0 otherwise. Here, we use a slightly 

different implementation from [53]. Indeed, using (2.14) we have 

„ p-i p p-i 

/ QdxUhdx(l))dx = ^ / (jJ-jdx4>)dx = ^ 'dKj&Vth, 


l—m 


l—rrij 


where the coefficients Dj are defined by dxUh = X]f=i Using 

the matrix M defined in (2.15), we get 


£>'■ = - ^ Mn,lUf, WO <l <p, Kj G H/i. 




We propose to apply this method in the context of our numerical experiments by 
selecting the range of modes where the spectral viscosity is applied according to the 
local smoothness of the solution: we set rrij = 1 when the solution is irregular and 
rrij = p — 1 otherwise. We apply the shock detection technique from m to test the 
local smoothness of the solution. 

Figure [6^ presents the convergence histories and solutions for the LLF numerical 
fliDc and third and fourth approximation orders. Results obtained with the above 
method and a viscosity coefficient e = h/{p + 1) are compared to that obtained 
without artificial dissipation. We observe the stabilization of the computations by 
the numerical viscosity. The amplitude of the oscillations are also lowered and are 
effectively damped in the uniform region where viscosity is applied only to the highest 
DOFs. The technique allows the computations to converge to steady-state solutions 
and to lower oscillations while keeping accuracy in smooth regions of the flow. 

7. Concluding remarks. Discrete shock profiles for a scalar conservation law 
with a convex flux discretized with a DG method have been analyzed. Using the 
Godunov numerical flux, we show existence of stationary profiles that are oscillating 
for polynomial degree p > I into one discretization cell only. The oscillations may 
vanish when the shock is located at an interface of the mesh. A linear stability 
analysis of the shock profiles show however that these latter solutions may be unstable. 
Gonsidering the inviscid Burgers equation, these profiles are constructed analytically 
for p < 3 and are shown to be parametrized by the shock strength and its relative 
position in the cell. 

The extension of this analysis to other numerical fluxes is also investigated. A 
theoretical analysis shows that oscillations propagate in neighboring cells but decay 
exponentially fast from the shock position for some class of monotone numerical fluxes. 
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Moreover, numerical experiments indicate that the shock profiles present strong simi¬ 
larities with the profiles obtained with the Godunov flux which may be considered as 
a relevant model for the analysis of the DG method. 

Finally, these results show that, when using a hierarchical functional basis, only 
the highest DOFs are responsible, at first order, of linear instability and propagation 
of oscillations in neighboring cells. As an application of this property, the spectral 
vanishing viscosity method is successfully used to stabilize computations and damp 
oscillations through a selective action on DOFs. 

One main contribution of this work is the analysis of oscillations and unstable 
character of the DG method for scalar conservation laws. These results may help to 
design specific stabilization techniques and future investigations will focus on these 
methods. 
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Appendix A. In this appendix, we prove the results of Theorem |4.1| which gives 
the solution of the numerical scheme in the cell where the shock is located as a function 
of the relative shock position, s 
first DOF = ulSc from (|3.8[) in Lemma 


c, given by (| 3.9[). We note that the solution for the 
We also give the following result 


3.4 


which holds for the Burgers equation and will be used to evaluate the local residuals 
associated to the equation for the second DOF: 


I f{uh)dx4’]dx = ^ y 


P P 


(A.l) 


1=0 


1=0 


where we have used the orthogonality of the function basis (2.91. According to the 
assumption of Theorem 4.1 the exact shock position is assumed to satisfy X- 1/2 < 
Xc < Xii 2 - From Theorems 3.1 and 3.2 it follows that the trace at the left interface 
satisfies = ul and a similar relation holds at the right interface Uy 2 — '^R- 

Further assuming (3.2), the numerical fluxes at the left and right interfaces of cell 
jc = 0 read h{uZ-^/ 2 > '^- 1 / 2 ) = /('“I 1 / 2 ) = and h{uZ/ 2 , ut/ 2 ) = fi'^t/ 2 ') = /(^r)- 

As a consequence, the equation for the first DOF always reduces to the trivial relation 


/(rl) = /(rr)- As suggested in Theorem 4.1 we adopt the following notatio n fo r 
the DOFs in the cell: Uq = ului, for all 0 < Z < p, where uq = Sc according to (3.8). 


Finally, we observe that according to (2.11), assumption (3.2) may be rewritten as 
follows 


^- 1/2 

UL 


=j2(-^yui > - 1 , 

1^0 


?/“ p 
Ul 


(A.2) 


A.l. solution for p 

DOF reads 


= 1. Using (A.l), the numerical scheme for the second 


— Mi(uQ + -g!") -I- /(ml) + /(m_r) — m| ^ —3^ ~ (^-3) 


whose solution reads ui = ±-\/3(l — s^). The condition (A.2 r) rea ds Sc —mi > —1 and 
restricts the solution mi > 0 over (f, 1), while the condition (A.2 o) reads Sc + Mi < 1 
and restricts the solution mi > 0 over (— 1, —^ ). The negative solution is the only one 
valid over (—1,1) which gives the result (|4.l[). 
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A.2. solution for p = 2. Equations for the second and third DOFs read 

- / f{uh)dx(t)odx + f{uL) + f{un) = 0 , 

J K 

- [ f{uh)d^4>ldx + f(uL) - f{uR) = 0, 

J K 


and give 


-Se-^-^ + 1 = 0 , 


-2mi(Sc+^W2) = 0. 


(A.4a) 
(A.4b) 


The second equation has a trivial solution ui = 0 which gives U 2 = ±\/5(l — s^). 
The validity of these solutions is imposed by (A.2) which reduces to —1 < s ± 
•^5(1 — s^) < 1 because the left and right traces are identical. Therefore, the so¬ 
lution U 2 = a/ 5(1 — Sg) is valid when s < — | an d th e solution U 2 = — s^) is 

valid when s > | which corr espo nd to solutions (4.2 3,,c). 


The second solution of (A.dj) reads U 2 = — |sc and then (A.4j,) leads to ui = 


Again, conditions (A.21 impose the 


± 1 / 3(1 — 9s2/4) which are valid only if |sc| < 
negative solution m = —•/3(1 — 95^/4). 

A.3. solution for p = 3. The numerical scheme for the second, third and fourth 
DOFs give 


„2 “2 

^^3 5 

M 3 

- 7 

= 0 , 

(A.5a) 

„ ^ 2 , 
-2 mi(Sc ± -M 2 j - 

18 

- 1^U2U3 

35 

= 0 , 

(A.5b) 

9 9 9 Uo G 

M 1 M 3 ± 1 

= 0 . 

(A.5c) 


We first observe that for a solution {ui,U 2 ,U 3 ) to (A.5), then (—ui,U 2 ,—U 3 ) is 
also solution. Therefore, we only look for solutions with mi < 0 and deduce the other 
ones by symmetry. Then, we also note that the choice U 2 = 0 imposes rti = 0 through 

I y I 2 2 

(A .5 d) and then (A.5 r,c) reduce to 1 —— ^ =0 and 1 —— ^ =0 which is possible 


only if M 3 = 0 and Sc = ±1 but is excluded by the s trict inequalities in (A.2). In the 
following, we thus consider U 2 7 ^ 0. Equation (A.5 d) induces M 3 = — ^(5Sc ± 2 m 2 ) 


and subtracting 7x(A.5i) and 3x(A.5:), multiplying by M 2 , one obtains 


-10(^ + Sc)ul + ^M 2 ± 6 scM 2 ± 4(1 - s^)m2 = 0, 
o oO 


whose solution for mi reads 


ul = 


3m2 


5(m 2 ± 3Sc) 
and one deduces the solution for M 3 from 

«3 = -5^(5s. 


+ 3ScM2 ± 2(1 — , 


2 m 2 


(A. 6 ) 


(A.7) 
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which gives 


2 49(5Sc + 2u2)^ , o , on 2^^ 

> = imu,(u, + 3s,) (35 + + 2(1 - * 0 )j^ 


Substituting the last result and (A. 6 ) into (A.5i), one obtains the equation for 
U 2 only: 


2(7Sc + U 2 ) (^iul + IScul + 14s^U 2 - + 7sc(l - = 0. 


(A. 8 ) 


Taking symmetries into account, the system (A.51 has at most height real solu¬ 


tions. The first root of equation (A. 8 ) U 2 = —7sc leads to the first solution (4.3i): 


ui = -U 3 = ±y V - 54s2, 


if |sc| < y ^ from which we again retain the solution ui < 0. Assumption (A.2 1 


requires m^_i /2 = ~7Sc > —1 — Sc and it^+ 1/2 = ~7sc < 1 — Sc which is satisfied if 


and only if — ^ < Sc < I- The three other roots of equation (A.5) read 


U 2 = 


U 2 = 


U 2 = 


1 

12 . 

1 

12 

1 

12 


- 7Sn - 


7(17_4- 6 ) , ^ 1/3 


A, 


1/3 


^ lillsl - 6) ^ 1/3 

vri/3 ^3 


a: 


- 7Sn - 


7(l-f73)(17sg-6) _ (1 + z73) ^i/3 


2(1 




(A.9a) 
(A.9b) 
(A.9c) 


where A 3 = 7 A 2 , A 3 = — 7 A 2 and A 2 has been defined in equation (4.4i). The last 
root (A.9:) is complex. From (4.4), it follows that the first root ( A. 9^ ) exists only 
when Aj = 7776s(! — 100085 ^ -|- 33593 ^ — 56 > 0 and A 2 > 0. The former condition 
has only two real roots Sc = ±si and is satisfied providing that |sc| > si defined by 


Si = 


18^2 L 


278- 


21/38411 -f (1373963 - 6687^15603)2/3 


2-1/3(1373963 - 6687\/lM)3)i/3 


while the latter condition has one root Sc = —\/6/17 and requires the additional 
conditions Sc < or Sc > Si. After substituting the former solution (A.9i) 

into (A. 6 ) and (A.7) one o btains the solution (4.3 1 ) for ui and U 3 when A 2 > 0. It 
may be checked that ( A.2 3 ) is always satisfied over the range [—1, —a/ 6/17] U [si, 1), 
while the condition (A.2r) is satisfied only over the range (—1, —-/6/17] U (1/6,1]. 
The solutions ui and M 3 with opposite signs cannot satisfy those conditions for any s 
value. 


Likewise, the second solution (A.9 d) hol ds w hen A 2 = —A 2 > 0, that is wh en 
-a/6/17 < Sc < Si and leads to solution ( |4.3p ). This solution satisfies ( |A.2| if 

-i/^<Sc < -1/6. 
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(a) Godunov 




(c) Osher 


Fig. 6.7. Eigenspectra of the linearized operator (5.3^ around steady-state solutions of the 
Burgers equation for different kind of numerical fluxes (GOD: Godunov, LLF, and OSH: Engquist- 
Osher) and 0 < Sc < 1, using the forward Euler (RKl) and explicit third-order Runge-Kutta (RK3) 
time integration schemes with 1 < p < 2, = 20, and A = l/(2p + 1) (open symbols) or X = 

0.1/(2p + 1) (full symbols). The unit circle 0(0,1) denotes the stability domain. 
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p=2,s^=7'10,GOD(^=10.217) - p=2, s,=4/5, GOD (ji=110.412E} 




(a) Sc = ^ (b) Sc = I 


Fig. 6.8. Real part of eigenvector associated to most amplified eigenvalues of the linearized 
operator \5.3]) around steady-state solutions of the Burgers equation for different kind of numerical 
fluxes (GOD: Godunov, LLF, OSH: Engquist-Osher), using the forward Euler (RKl) and explicit 
third-order Runge-Kutta (RK3) time integration schemes with p = 2, N = 20, and A = 






(a) 


(b) 


Fig. 6.9. Oonvergence histories (top) and final solutions (bottom) without (left) and with the 
spectral viscosity method (right) of the Burgers equation using the LLF flux with p = 2 and p = 3. 

















































